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We study numerically the dynamical system of a two- 
electron atom with the Darwin interaction as a model to in- 
vestigate scale-dependent effects of the relativistic action-at- 
a-distance electrodynamics. This dynamical system consists 
of a small perturbation of the Coulomb dynamics for energies 
in the atomic range. The key properties of the Coulomb dy- 
namics are: (i) a peculiar mixed-type phase space with sparse 
families of stable non-ionizing orbits and (ii) scale-invariance 
symmetry, with all orbits defined by an arbitrary scale pa- 
rameter. The combination of this peculiar chaotic dynam- 
ics ((i) and (ii)), with the scale-dependent relativistic correc- 
tions (Darwin interaction) generates the phenomenon of scale- 
dependent stability: We find numerical evidence that stable 
non-ionizing orbits can exist only for a discrete set of resonant 
energies. The Fourier transform of these non-ionizing orbits 
is a set of sharp frequencies. The energies and sharp frequen- 
cies of the non-ionizing orbits we study are in the quantum 
atomic range. 

PACS numbers: 05.45+b, 31.15.Ct, 03.20.+i, 05.45.Pq 



The Coulomb dynamical system of the helium atom 
is a very peculiar chaotic system that exhibits Arnold 
diffusion Q, and with a typical trajectory having an in- 
finity of possible time-asymptotic final states. For exam- 
ple, almost all negative-energy trajectories of Coulom- 
bian helium display the generic phenomenon of ioniza- 
tion, namely, the ejection of one electron 0. Several 
nonlinear dynamical systems share this property of hav- 
ing more than one time-asymptotic final state, with the 
respective basins for each outcome having a complicated 
structure in initial condition space |^,^|. The numerical 
work on this paper is based on stable Coulombian orbits 
of a two-electron atom that do not ionize for several mil- 
lions of turns of one electron around the nucleus. It is 
a property of the Coulomb dynamics of a two-electron 
atom that most initial conditions with a negative energy 
ionize very quickly in about 20 turns [0 . Then there are 
the very special initial conditions that do not ionize due 
to a precise phase balance between the two electrons. 
These rare non-ionizing orbits are defined very sharply 
in phase space and were first studied in reference Q for 
plane orbits. Here we also develop a numerical procedure 
to search for non-ionizing orbits among a large number 
of possible tridimensional initial conditions. 

The Coulomb Hamiltonian exhibits the scale invari- 



ance degeneracy: if we scale time and space as t — » Tt , 
r — > Lr, for T 2 /L 3 — 1, the equations of motion are left 
invariant. For this reason, the behavior of the Coulomb 
dynamics is the same in all scales, a degeneracy which is 
broken by the relativistic effects of electrodynamics. The 
phenomenon of breaking the scale invariance in electro- 
dynamics was explored analytically in |^,^| for the Dar- 
win interaction, which is the low-velocity approximation 
to the Wheeler-Feynman action-at-a-distance electrody- 
namics 0. It was found in [js],^) that a simple resonant 
normal form approximation theory predicts a discrete set 
of quantized scales very close to the quantum atomic 
energies. Using these preliminary findings as guide, we 
present a numerical investigation of the stability of non- 
ionizing orbits for the Darwin dynamics and its depen- 
dence on the energy scale. It turns out that for energies 
of atomic interest, the Darwin equations of motion ap- 
proximate the Coulomb equations plus a perturbation of 
size /3 2 , with (3 ~ 10~ 2 . Therefore, non-ionizing stable 
orbits of the Darwin dynamics should exist in the neigh- 
borhood of non-ionizing stable Coulombian orbits if the 
perturbation does not force ionization. For these, our nu- 
merical results with the Darwin dynamics indicate that 
the non-ionizing property plus stabilty require sharply 
defined discrete energies. 

The Darwin interaction is not exactly a Lorentz invari- 
ant interaction [p|-^0| , so we study it as an approximation 
to the relativistic action-at-a-distance electrodynamics, 
for the sake of including the present approach into an un- 
derlying physical theory. Maxwell's theory would seem to 
be the natural candidate for the comprehensive physical 
theory, but it lacks time-reversibility and dipolar dissipa- 
tion would forbid the orbits studied in this paper. There 
is also the choice of other more recent Lorentz-invariant 
Lagrangian fi"l|| systems and constrained Hamiltonian 
dynamical systems Jl2|-p[, whose exact forms are ac- 
tually more amenable to numerical treatment than the 
Wheeler-Feynman electrodynamics, but we shall not con- 
sider them here. The interested reader should consult ref- 
erence [^2| , where a covariant approximation to Wheeler- 
Feynman electrodynamics is attempted by the two-body 
Todorov equation of constraint dynamics. 

This paper is organized as follows: in section I we re- 
view the state of the art of the time-reversible action- 
at-a-distance electrodynamics, and if the reader wants to 
skip this part the rest of the paper makes full sense as a 
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nonlinear dynamics study, except for the discussion at the 
end. In section II we describe the numerical calculations 
with the Coulomb limit of the Darwin interaction, find 
some non-ionizing orbits and their Fourier transforms. In 
section III we include the scale dependent Darwin terms 
and investigate the possibility of stable non-ionizing or- 
bits. In section IV we put the conclusions and discussion. 



I. ACTION- AT- A-DISTANCE 
ELECTRODYNAMICS 

The Wheeler- Feynman Jlq] electrodynamics developed 
from the Schwarzschild-Tetrode-Fokker [jl6| direct- 
interaction functional. Equations of motion are derived 
from Hamilton's principle for the action integral 
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where the four- vector Xi(st) represents the four-position 
of particle i parametrized by arc-length Si , double bars 
indicate quadri- vector modulus \\xi — Xj\\ 2 = (xi — Xj) ■ 
(xi — Xj) and the dot indicates the usual Minkowski rel- 
ativistic scalar product of four-vectors, (integration is to 
be carried over the hole particle trajectories, at least for- 
mally). The above action integral describes an inter- 
action at the advanced and retarded light-cones with 
an electromagnetic potential given by half the sum of 
the advanced and retarded Lienard-Wierchert potentials 
JtJ. Wheeler and Feynman showed that electromagnetic 
phenomena can be described by this direct action-at-a- 
distance theory in complete agreement with Maxwell's 
theory as far as the classical experimental consequences 
|[^,[L7|. This direct-interaction formulation of electro- 
dynamics was developed to avoid the complications of 
divergent self-interaction, as there is no self-interaction 
in this theory, and also to eliminate the infinite number 
of field degrees of freedom of Maxwell's theory It 
was a great inspiration of Wheeler and Feynman in 1945, 
that followed a lead of Tetrode |Q and showed that with 
the extra hypothesis that the electron interacts with a 
completely absorbing universe, the advanced response of 
this universe to the electron's retarded field arrives at 
the present time of the electron and is equivalent to the 
local instantaneous self-interaction of the Lorentz-Dirac 
theory . The action-at-a-distance theory is also sym- 
metric under time reversal, as the Fokker action includes 
both advanced and retarded interactions. Dissipation in 
this time-reversible theory becomes a matter of statisti- 
cal mechanics of absorption p0| . The area of Wheeler- 
Feynman electrodynamics has been progressing slowly 
but steadily since 1945: Quantization was achieved by 
use of the Feynman path integral technique and the ef- 
fect of spontaneous emission was successfully described in 
terms of interaction with the future absorber, in agree- 
ment with quantum electrodynamics |21j. It was also 



shown that it is possible to avoid the usual divergen- 
cies associated with quantum electrodynamics by use of 
proper cosmological boundary conditions pl| . As far as 
understanding of the dynamics governed by the equa- 
tions of motion, the state of the art is as follows: The 
exact circular orbit solution to the attractive two-body 
problem was proposed in 1946 |22| and rediscovered by 
Schild in 1962 p3[ . The 1-dimensional symmetric two- 
electron scattering is a special case where the equations 
of motion simplify a lot and it has been studied by many 
authors, both analytically and numerically [p4]-p6[. In 
this very special case the initial value functional problem 
surprisingly requires much less than an arbitrary initial 
function to determine a solution manifold with the ex- 
tra condition of bounded manifold for all times. It was 
shown that the solution is uniquely determined by the 
intcrclcctronic distance at the turning point if this dis- 
tance is large enough (this minimum distance curiously 
evaluates to 0.49 Bohr radii by the action-at-a-distance 
theory pjj, much larger than about one classical elec- 
tronic radius that one would naively guess). As a result 
of this theorem, there is a single continuous parameter 
(the positive energy) describing the unique non-runaway 
symmetric orbit at that given positive energy. 

The Noether's four-constant of motion derived from 
the Fokker Lagrangian involves an integral over the past 
history [[7| JT5|j2^ | . For example in the case of a hydrogen 
atom this four-momentum constant evaluates to 

P x = m p Xp + eA x (x p ) + m e x x - eA x (x e ) 

26 / d/Tp I dT e S ( \\x e -EpW J X e )x e Xp 
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+2e 2 / dr p I drj (\\x e - x p \\ 2 ) (x p - x e )x e x p , 
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where 5 represents the derivative of the delta function 
0. Notice that because of this delta function, only fi- 
nite portions of the trajectory are involved: actually an 
extent of length t ~ Iryijc approximately. This non- 
local constant will behave very differently from the local 
Coulombian energy, that is known to confine orbits of a 
negative energy within a maximum separation distance. 
In the case where the particles acquire a large separa- 
tion (unbound state), the hole past history is involved 
(i ~ 2ri2/c ~ oo) in the determination of the non-local 
energy constant. 

As regards the mathematical structure of the equa- 
tions of motion, for the case of a two-electron atom the 
acceleration of electron 1 is given by P7J 
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where — e and m are the electronic charge and mass, 7 = 
1/^/(1 _ llWy^and E and B 

are the total electric and 
magnetic fields produced by electron 2 and the nucleus. 
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In the action-at-a-distance theory these fields are given 
by the average of the retarded and advanced Lienard- 
Wiechert fields, calculated with the instantaneous posi- 
tion of the stationary nucleus and the retarded and ad- 
vanced positions of electron 2 at the times ti = , which 
is defined by the implicit condition 



\r 2 (t T )-r 1 (t 1 )\=Tc(t T ^t 1 ), 



where the minus and plus signs are the conditions for the 
retarded and advanced times respectively. The partial 
electric fields of electron 2 acting on electron 1 at time 
h are pi 
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where R T n T = r 2 (t T ) - n(ii) , f3 2 = v 2 /c,(-e),j 2 = 
(1— /3|) _1//2 and c is the speed of light. The advanced field 
E+ is obtained from the above expression by replacing t_ 
by t + and c by — c. The partial magnetic fields of electron 
2 are 



±n T x E T , 



where the ± is to ensure an outgoing Poynting vector 
(cE x B) for the retarded fields and an incoming Poynt- 
ing vector for the advanced fields. The total electric 
field in equation (Q) must include also the instantaneous 
Coulomb electric field of the stationary nucleus. 

Equation (|l|) can suggest a paradox about causality, 
as the force depends on the future of particle 2. In the 
following, and to finish this introduction, we show that 
equation (Q), when written properly, becomes a func- 
tional differential equation with delayed argument only, 
as first observed in J29] . To outline the essentials of the 
explanation, let us first ignore the field of the nucleus 
and take the nonrelativistic limit of (0) (vi =0). In 
this approximation the electric field E entering in equa- 
tion (l|) evaluates to E = 0.5E+(xi, x 2 +,v 2 +, a 2 +) + 
0.5E- (xi,x 2 -,v 2 -,a 2 -). Then we note that one can use 
equation ([!]) as an equation of motion for particle 2 , by 
solving the rearranged form of (|l]) , 

eE + (xi,x 2+1 v 2+ ,a 2+ ) = 

—2mai(t) - eE-(xi,X2-,V2-,a2-), 

(2) 

for the most advanced acceleration of particle 2, a 2+ = 
a 2{t + — )■ In the above form it is clear that the right 
hand side involves only functions evaluated at times prior 
to the most advanced time, defined by s = t ^ — — , and 
no further advanced information is necessary, eliminat- 
ing the ghost of dependence on the future. In the same 



way, the causal equation of motion of particle 1 is to be 
produced from the equation for particle 2 by solving for 
the most advanced acceleration of particle 1. For the 
special case of 1-dimensional motion of two electrons, 
E+ = E + (xi,x 2+ ,v 2+ ) depends only on the advanced 
velocity, and (||) can easily be solved for this advanced 
velocity as a function of the past history. In the 3- 
dimensional case there is an extra complexity, as the ac- 
celeration appears in the Lienard-Wicchert partial field 
E+ in the form n\ 2 x (n\ 2 x a 2 +)/ri 2 . The bad news 
is that the component of the acceleration along the ad- 
vanced normal can not be solved for from the value of the 
double- vector-product only. Because of this degeneracy, 
equation (§) is an algebraic-differential equation, and the 
null direction of the left hand side of @) is a constraint 
to be satisfied by the right hand side (the scalar product 
with ri\ 2 must vanish). The numerically correct way to 
integrate this type of equation is by use of the modern in- 
tegrators for algebraic-differential equations like DASSL 
fiOf adapted for retarded equations (which has never been 
done yet) or by dealing directly with the algebraic con- 
straint |3l]]. According to the standard classification of 
G. A. Kamenskii p2]| , equation ^ belongs to the class 
of differential-difference equations of neutral type. Even 
though more complex, the motion is still causally deter- 
mined by the past trajectory, as we wanted to demon- 
strate, the price being an algebraic neutral delay equa- 
tion. 

As far as initial conditions go, the general theory on 
delay equations |3^] tells us that we need to provide an 
initial C 2 function describing the position of particle 2 
from s — — t — ^= up to the initial instant s = 0. 

c c L 

The information on particle 1 needed is also to be pro- 
vided over twice the retardation lag seen by particle 1. 
This is a short piece of trajectory for bound nonrela- 
tivistic atomic orbits, but for a ionized state or a run- 
away orbit this can be the whole past history! Unless 
further simplifications or conditions are added, this is 
the generic problem at hand. The 3-dimensional cases 
of atomic interest (e.g. helium) have never been stud- 
ied, and they are more complex than the 1-d scattering 
because one can have negative energy bound states for 
example. Most relevant for physics is the question of the 
conditions for the existence of a bounded manifold solu- 
tion, which still needs to be understood in the general 
case (it would be very curious if they turned out to be 
a discrete set of negative energies). The only existing 
analytical result in the 3-dimensional case is the linear 
stability of the Schonberg-Schild circular orbits ^7j, re- 
sulting in an infinite number of unstable solutions to the 
characteristic equation. The numerical treatment ofthe 
exact neutral equations displays instabilities and is gen- 
erally difficult. In the following we resort to the Darwin 
approximation not as much as a mathematical approxi- 
mation to the action-at-a-distance electrodynamics, but 
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as a physical approximation of Lorentz-invariant dynam- 
ics in the atomic (shallow) energy range. 



II. NUMERICAL CALCULATIONS FOR THE 
COULOMB DYNAMICS 

To introduce our numerical calculations, we start from 
the scale-invariant Coulomb limit of the Tetrode- Fokker- 
Wheeler-Feynman interaction: Let — e and m be the elec- 
tronic charge and mass respectively and Ze the nuclear 
charge of our two-electron atom, which in this work is as- 
sumed to have an infinite mass. All our numerical work 
uses a scaling which exploits the scale invariance of the 
Coulomb dynamics: Given a negative energy, there is a 
unique circular orbit at that energy with frequency u) a 
and radius R related by e 2 /(mui 2 R 3 ) = 1/(Z - |) = 
C(Z). We scale distance, momentum, time and energy as 
x — > Rx, p — > mu> Rp , L0 o dt — ► dr and E — > mLu 2 R 2 H, 
respectively. In these scaled units, the Coulomb dynam- 
ics of the two-electron atom is described by the scaled 
Hamiltonian 



H=i(lft| 2 + |ft| 



<(Z){-----}, (3) 
m n r 2 



where n = \xx\, r 2 = |s 2 |, r*i2 = |xi — x%\ (single bars 
represent euclidean modulus) and (3 = lo R/c. For a 
generic non-circular orbit, (3 plays the role of a scale pa- 
rameter, and we recover the value of the energy in ergs 
through E = mc 2 f3 2 H. Notice that (3 does not appear 
in the scaled Hamiltonian, which is the scale invariance 
property. From the scaled frequency w and scaled angu- 
lar momentum I we can recover the actual values in CGS 
units by the formulas 



mc 2 C(Z)/3 3 
e 2 /c 



I = 



e 2 /c I 

az)T 



(4) 



The only other analytic constant of the Coulomb dynam- 
ics, besides the energy (g) is the total angular momen- 
tum, and this dynamics in chaotic and displays Arnold 
diffusion, as proved in jjj for a similar three-body system. 

The numerical calculations were performed using a 
9th-order Runge-Kutta embedded integrator pair [f33"| . 
We chose the embedded error per step to be 10~ 14 , and 
after ten million time units of integration the percent- 
age changes in energy and total angular momentum were 
less than 10~ 6 . As a numerical precaution we performed 
the numerical calculations using the double Kustanhcimo 
coordinate transformation to regularize single collisions 
with the nucleus |j4j . As these alone are not enough for 
faithful integration, we checked that there was never a 
triple collision, as the minimum inter-electronic distance 
was about 0.3 units while the minimum distance to the 
nucleus was 0.01 units for all the orbits considered in this 
work. We also checked that along stable non-ionizing or- 
bits we can integrate forward up to fifty thousand time 



units, reverse the integration, go backwards another fifty 
thousand units and recover the initial condition with a 
percentile error of 10 -5 . For longer times this precision of 
back and forth integration degenerates rapidly, which is 
due to the combined effect of numerical truncation and 
stochasticity. The question of how far in time the nu- 
merical trajectories approximate shadowing trajectories 
in the present system is far from trivial p5| |, but we as- 
sume it to be a time at least of the order of these one 
hundred thousand units. (Energy conservation of one 
part in a million is achieved for much longer times, even 
one billion time units). 

The study of orbits of a two-electron atom was greatly 
stimulated by the recent interest in semiclassical quan- 
tization, and these studies discovered two types of sta- 
ble zero-angular-momentum periodic orbits for helium 
(Z = 2): the Langmuir orbit and the frozen-planet or- 
bit |36| , |37|| . A detailed study of the non-ionizing orbits 
of Coulombian helium was initiated in reference Q for 
plane orbits, and we describe some of their results be- 
low. There are basically two types of non-ionizing or- 
bits: Symmetric if r\ = r 2 for all times and asymmet- 
ric if n 7^ r 2 generically. Symmetric orbits are pro- 
duced by symmetric initial conditions like for example 
Xl (0) = -x 2 (0) and «i(0) = -v 2 (0) or ^(0) = -x 2 (0) 
and «i(0) = u 2 (0) with Xl (0) ■ ui(0) = §|] Because (§) 
is symmetric under particle exchange, these orbits satisfy 
ri = r 2 at all times, and therefore cannot ionize if H < 
(both electrons would have to ionize at the same time, 
which is impossible at negative energies). For example 
the double-elliptical orbits (two equal ellipses symmetri- 
cally displaced along the x-axis) discussed in Q are in 
this class. Double-elliptical orbits are known to be un- 
stable |5|,|| and we find that they ionize in about one 
hundred turns because of the numerical truncation error. 
Most symmetric plane orbits are very unstable to asym- 
metric perturbations, with the exception of the Langmuir 
orbit for a small range of Z values around Z = 2 |3(| . 

The simplest way to produce an asymmetric non- 
ionizing plane orbit is from the initial condition X\ = 
(ri,0,0) , x x = (0,1-1^4/7,0) , x 2 = (-1.0,0,0) , 
%2 = (0, — v4~A> 0), as suggested in 0. In Figure 1 we 
show the electronic trajectories for the first three hundred 
scaled time units along a two-dimensional non-ionizing 
orbit of Ca+ ls (Z = 20) with r x = 1.4 and vx = 1.28442 
in the above defined condition. We used a numerical re- 
fining procedure to finely adjust wias to maximize the 
non-ionizing time and this condition of Figure 1 does not 
ionize for one million time units. The orbit survives that 
far only for a very sharp band of values of vx , other neigh- 
boring values producing quick ionization. This orbit was 
named double- ring torus in ||. The other possible type 
of non-ionizing orbit resulting from the above initial con- 
dition, depending on (rx, Vx), is what was named braiding 
torus in reference 0] , with both electrons orbiting within 
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the same region. A search over (n, v\) was conducted in 
||| , and it was found that most values of (ri , «i ) produce 
quick ionization except for a zero-measure set of (ri,i>i) 
values where braiding tori or double ring orbits are found. 
This suggests the general result that non-ionizing orbits 
are rare in phase space. 

To search for general tridimensional non-ionizing orbits 
in phase space, it is convenient to introduce canonical 
coordinates Xd and x c 

Pd = (pi -p2)/V2, x d = (xi - x 2 )/V2 

p c = (pi +p 2 )/V2, x c ee (xi + x 2 )/V2. (5) 

Initial conditions with x c = p c = describe double- 
elliptical orbits (and circular as a special case). To gen- 
erate an elliptical initial condition, we exploit the scale 
invariance and set the energy to minus one. It is easy to 
check that elliptical orbits of the Hamiltonian (||) with an 
energy of minus one must have a total angular momen- 
tum of magnitude ranging from zero to two. To exploit 
the rotational invariance of (||), we can choose the plane 
defined at x c — p c = by the angular momentum L = 
Xd x pd + x c x pc = Xd x pd to be the xy plane. On 
this xy plane a single number < \xa X pd| < 2 (the 
angular momentum), determines completely the ellipti- 
cal orbit. The next step in producing a generic orbit 
is to add all possible perturbations along x c and p c to 
the chosen elliptical orbit. These are six directions and 
once we are looking for bound oscillatory orbits, we can 
choose z c = 0, once z c has to cross the xy plane at some 
point. These are five numbers to vary and plus the angu- 
lar momentum of the elliptical orbit it totals six parame- 
ters. Our numerical search procedure consists in varying 
these six parameters over a fine grid, integrating every 
single initial condition until the distance from one elec- 
tron to the nucleus is greater than twenty units, which is 
our ionization criterion. This criterion fails if the orbit 
has a very low angular momentum because these can go 
far away from the nucleus and come back, and therefore 
our search possibly misses low-angular-momentum non- 
ionizing orbits. As the majority of the initial conditions 
ionize very quickly, this search procedure is reasonably 
fast. We first perform a coarse search for ionization times 
above one thousand units and then refine in the neighbor- 
hood of each surviving condition to get conditions that 
do not ionize after one million time units. 

Using the above numerical search procedure we found 
the tridimensional non-ionizing initial condition of Figure 
2 for helium, a tridimensional double-ring orbit generated 
by the initial condition 

x x = (1.2812617,0.0147169,0.0) 

x 2 = (-1.5511484,0.0147169,0.0) 

pi = (-0.0194868,0.4398889,0.1094930) 

p 2 = (-0.0194868,-0.7972467,0.1094930), 



which does not ionize before ten million turns. (After 
the search and refinement, we scaled this orbit's energy 
to minus one, for later convenience). We also found the 
non-ionizing orbit orbit of Figure 3 for H-minus [Z = 1), 
a tridimensional orbit generated by the condition 

xi = (1.9776507,-0.3411364,0.0) 
x 2 = (-1.2288121,-0.3411364,0.0) 
pi = (0.0421302,0.5057782,0.2810539) 
p 2 = (0.0421302,-0.4132970,0.2810539), 

which does not ionize before one million turns (Coulom- 
bian energy of this condition is also minus one). This 
last orbit is fragile and numerically harder to find: as the 
first electron has an orbit very close to the positive Z = 1 
charge, there remains only a dipole field to bind the sec- 
ond electron. As the outer electron is much slower in the 
scaled units, we had to plot the first 10000 time units of 
evolution to display the generic features of the trajectory. 
Non-ionizing orbits of H~ are very rare in phase space, 
which is reminiscent of the quantum counterpart, as the 
H~ ion is known to have only one quantum bound state 
at E ~ — 0.55mc 2 a 2 , very close to the ionization thresh- 
old (-0.5mc 2 a 2 ) ||. 

One remarkable fact about these non-ionizing orbits is 
that they all have a very sharp Fourier transform. This 
property makes them approximately quasi-periodic or- 
bits. For example in Figure 4 we plot the fast Fourier 
transform of the orbit of Figure 2, performed using 2 16 
points. (It seems that there are at least two basic fre- 
quencies in the resonance structure of Figure 4). Even 
though these orbits look like quasi-periodic tori, there 
seems to be a thin stochastic tube surrounding each or- 
bit, as evidenced by a small positive maximum Lyapunov 
exponent. We calculated numerically this maximum Lya- 
punov exponent by doubling the integration times up to 
T = 10 7 and found that the exponent initially decreases 
but then saturates to a value of about 0.001 for the or- 
bits of Figures 1, 2 and 3. The gravitational three-body 
problem has recently been proved to display Arnold dif- 
fusion [Q , and this numeriacally calculated positive Lya- 
punov exponent suggests that the same is true for the 
two-electron Coulombian atom. 

III. NUMERICAL CALCULATIONS FOR THE 
DARWIN DYNAMICS 

The numerical integrations in this section are per- 
formed using the Darwin approximation. The Darwin 
equations of motion are a 1 perturbation of the Coulomb 
dynamics, of size 1 ~ 10 -4 for atomic energies. In the 
scaled units of section II the Darwin Hamiltonian is the 
following (3 2 perturbation of Hamiltonian (||) 

Hz, = iflpii 2 + w 2 \ 2 ) + az){— ----} 

2 no 7*i r 2 
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r [pi -P2 + («i2 •pi)(ni2 • pb)] 

^12 

-|W + I#A (6) 

where fii2 = {x\ — 2?2)/?"i2' The second line represents 
the Biot-Savart magnetic interaction plus the first rela- 
tivistic correction to the static electric field and the last 
line describes the relativistic mass correction. Notice that 
these are both proportional to the small parameter /3 2 , 
which makes them a small scale-dependent perturbation 
on the scale invariant Coulomb Hamiltonian (first line). 
It is possible to regularize the Darwin equations with the 
same double- Kustanheimo transformation p4| ■ on ly that 
here one needs to define the regularized time using the 
higher powers dt — r\r\ds, instead of the lower powers 
dt = ri^ds used to regularize the Coulomb equations 

The main question we address numerically in this sec- 
tion is the dependence of the stability of a non-ionizing 
orbit with the energy scale of the orbit . Here we use 
the word stability to mean ionization-stability: We call 
an initial condition ionization-stable if any small pertur- 
bation of it produces another non- ionizing orbit. The 
scale-dependent Darwin terms (of size f3 2 ) produce sig- 
nificant deviations from the Coulomb dynamics only in 
a time-scale of order 1//3 2 , which we find numerically to 
be the typical time for a non-ionizing Coulombian initial 
condition to ionize along the Darwin vector field. This 
poses a numerical difficulty if (3 is too small because one 
has to integrate the orbit for very long times to inves- 
tigate the stability. It turns out that ionization-stable 
orbits can be found at larger values of j3 for larger val- 
ues of Z . Here the dynamical stability mechanism is 
reminiscent of quantum atomic physics, where the values 
of j3 vary with the nuclear charge as j3 ~ Z/137. Large 
values of Z facilitate the numerical procedure and in the 
following we present the numerical investigation of the 
stability of non-ionizing orbits starting from the large Z 
case. 

Let us start with the Z — 20 calcium ion two-electron 
system along the non-ionizing orbit of Figure 1 by fixing 
T\ = 1.4 and v\ = 1.28442. in the condition defined in 
section II. To test the stability of the orbit at each value 
of (3 we add a random perturbation of average size /3 2 
to the initial condition and integrate the Darwin dynam- 
ics until either we find ionization or the time of integra- 
tion is greater than 10 7 time units We repeat this for at 
least twelve randomly chosen perturbations (because of 
the twelve degrees of freedom) and the minimum time to 
ionization is plotted in Figure 5 as a function of (3. It 
can be seen that only for a narrow set of values around 
[3 ~ 0.037 this minimum time to ionization was greater 
than 10 6 . or the other values it decreases rapidly to a 
value of about 10 3 . One could argue that for the other 
values of (3 the non-ionizing initial condition has shifted 



away from the v\ = 1.28442 initial condition and this 
being the reason that our orbit ionized. To test this, 
we fixed (3 at a "bad " value for example (3 = 0.02 and 
varied the plane initial condition in the neighborhood of 
this condition of Figure 1 . We found that the minimum 
time to ionization was always about 10 3 (also the max- 
imum time before ionization was about 10 3 ). We also 
searched in a bigger neighborhood, of size proportional 
to (3. This suggests the interpretation that for the special 
resonant value of (3 = 0.037 the net diffusive effect of the 
scale-dependent term vanishes, allowing a non-ionizing 
perturbed manifold. In order to have a direct interpre- 
tation (in atomic units) of the scale parameter j3, it is 
convenient to scale to minus one the energy of the ini- 
tial condition of Figure 1 (by exployting the Coulom- 
bian scale invariance). After this, the energy of the orbit 
in ergs evaluates to E = mc 2 f3 2 H = —mc 2 f3 2 , and for 
(3 = 0.037 this is approximately —24.59 atomic units. 
The total angular momentum of this orbit is l z = 7.94?i. 
This orbit's energy is above the ionization continuum of 
the ion, E = -mc 2 a 2 Z 2 /2 — —200 atomic units, but it 
is still in the quantum range. It serves nevertheless to 
demonstrate that this dynamical system might exhibit 
non-ionizing stable orbits only at very sharply defined 
energy values. 

For the orbits of Figures 2 and 3, the above procedure 
becomes prohibitively slow, as the value of (3 are much 
smaller and one must integrate for very long times, much 
beyond the estimated shadowing time. To partially over- 
come this we used a larger amplitude random perturba- 
tion (of average size 20/3 2 ), to produce faster ionization. 
The drawback with this is that the minimum ionization 
time does not show pronounced peaks, only the average 
ionization time still showing a signature of scale depen- 
dence. In Figure 6 we show this average time for the orbit 
of Figure 3. This property of sharply defined energies can 
possibly be found for the lower-lying energies below the 
ionization threshold as well. These orbits would involve 
configurations where the electrons come very close to the 
nucleus and acquire a large velocity. Even though our 
integrator is regularized, the correct physical electronic 
repulsion is greatly amplified when one electron has a rel- 
ativistic velocity and the Darwin approximation can not 
describe the physics then. Actually, it is known that the 
Darwin interaction can produce unphysical effects when 
pushed to relativistic energies |4l). We therefore do not 
expect to find these low-lying atomic energy scales with 
the present Darwin approximation and shall be contempt 
with these interesting result already. 

For the same reason given above, we do not study here 
the frozen-planet periodic orbit (the two electrons per- 
forming one-dimensional periodic motion on the same 
side of the nucleus, with the inner electron rebounding 
from the origin, an artifact of regularization) . The main 
problem being the failure of the Darwin approximation, 
as the inner particle goes to the speed of light ElJ . The 
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correct relativistic dynamics can actually produce a new 
physical inner turning point very close to the origin but 
not at the origin as the regularized motion, and we dis- 
cuss elsewhere Q. 

Last, we consider the non-ionizing symmetric periodic 
orbit called the Langmuir orbit, where the two elec- 
trons perform symmetric bending motion shaped approx- 
imately like a semi-circle |59|. For the Coulomb two- 
electron atom with Z = 2 this orbit was found to have 
a zero maximum Lyapunov exponent [ |36| . The orbit 
is therefore neutrally stable, which is the best one can 
expect from a periodic orbit of a Hamiltonian vector 
field. (Absolute stability violates the symplectic sym- 
metry, which says that to every stability exponent A one 
should have a 1/A exponent). It is a simple matter to 
obtain the Langmuir-like orbit for the Darwin Hamil- 
tonian at any given value of ft: all it takes is a little 
adjusting in the neighborhood of the Coulombian Lang- 
muir condition. We attempted to investigate numerically 
any scale-dependent diffusion away from this Darwin- 
Langmuir condition for ft in the atomic range, but again 
the numerics is prohibitively slow at the time of writing 
this work. 

IV. CONCLUSIONS AND DISCUSSION 

The simplified dynamical mechanism behind resonant 
non-ionization seems to go intuitively as follows: The pe- 
culiar scale-invariant Coulomb dynamics determines the 
non-ionizing orbits within narrow "stochastic tubes ". 
The next step is the action of the small scale-dependent 
relativistic corrections that produce a slow diffusion of 
the orbit out of the thin tube in a time of the order of 
1/ ft 2 - After this, quick ionization follows. Only at very 
special resonant values of ft the relativistic terms leave 
the orbit within the tube, a resonant effect that depends 
on ft, fixing the energy scale. In the literature, the es- 
cape to infinity from simpler to understand two-degree- 
of-freedom systems has been attributed to cantori, which, 
as is well known, can trap chaotic orbits near regular re- 
gions for extremely long times In the present larger 
dimensional case it appears that resonances are also con- 
trolling the escape to infinity of one electron by the ex- 
istence of extra resonant constants of motion |||6| . This 
seems to be in agreement with the numerical results of 
very sharp peaks for the minimum ionization time. We 
have tried to concentrate on the physics described by this 
combination of chaotic dynamics on a two-electron atom 
with inclusion of relativistic correction, while discussing 
this highly nontrivial result of nonlinear dynamics. 

In references ||[| we noticed that a simple resonant 
normal form criterion gives a surprisingly good predic- 
tion for the discrete atomic energy levels of helium. The 
resonant structure was calculated using the Darwin in- 
teraction (ph, which is the low- velocity approximation 



to both Maxwell's and Wheeler-Feynman's elec- 
trodynamics. As we saw in section II, the Coulombian 
non-ionizing orbits are far from circular, and these or- 
bits would radiate even in dipole according to the time- 
irreversible Maxwell's electrodynamics (circular orbits 
radiate only in quadrupole but are linearly unstable). 
It becomes then clear that the heuristic results of 
can only have a physical meaning in the context of a 
time-reversible theory (as the action-at-a-distance elec- 
trodynamics for example). 

The combination of chaotic dynamics with relativis- 
tic invariance has never been explored numerically, and 
most known Lorentz-invariant dynamical systems are for 
one particle and possess trivially integrable dynamics. 
The situation gets unexpectedly much more complicated 
for more than one particle (apart from the trivial non- 
interacting many-particle system): Due to the famous 
no-interaction theorem K| , the relativistic description 
of two directly interacting particles is impossible within 
the Hamiltonian formalism and its set of ten canonical 
generators for the Poincare group JlJ] . Description of in- 
teracting particles is possible only in the context of con- 
straint dynamics, with eleven canonical generators and 
with the Dirac bracket replacing the Poisson bracket. For 
example the relativistic action-at-a-distance equations 
for two interacting electrons are non-local and possess 
only infinite-dimensional constrained Hamiltonian rep- 
resentations The interested reader should con- 
sult some recently found two-body direct-interaction rel- 
ativistic Lagrangian dynamical systems as well as the 
constraint-dynamics direct-interaction models recently 
used in chromodynamics and two-body Dirac equations 
p2| p"4| . The nonlinear dynamics of these models could 
display interesting and so far inexplored dynamical be- 
haviour. 

It would be natural to wonder if one can find an anal- 
ogous scale-dependent dynamics for a dynamical system 
describing the hydrogen atom, apparently the simplest 
example of Lorentz-invariant two-body relativistic dy- 
namics of atomic interest. It turns out that hydrogen 
is not simpler than helium at all, but it appears to us 
that there is an essential difference which has actually 
made the interesting dynamics of a two-electron atom 
amenable to study already within the Darwin approxi- 
mation: In a two-electron atom orbits with a negative 
energy can ionize, while in hydrogen this might be pos- 
sible only if one includes all orders of the relativistic 
action-at-a-distance interaction. (As we saw in section 
I, the "Noether's energy constant" involves a segment 
of the past trajectory, and a negative value does not for- 
bid ionization). Ionization with a negative energy would 
be impossible for hydrogen within the Darwin approxi- 
mation (unless the electron goes to the speed of light). 
This is indication that in hydrogen the essential physics 
described by the action-at-a-distance electrodynamics is 
of non-perturbative character. The paradoxical result of 
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the infinite linear instability of circular orbits in atomic 
hydrogen p?j] is another warning of this non-perturbative 
dynamics. 
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FIG. 1. Non-ionizing double-ring orbit for Ca + ion 
(Z = 20), obtained from the initial condition x\ = (n, 0, 0) , 
xi = (0,ui v /4/7,0) , x 2 = (-1.0,0,0) , ±2 = (0,-^/477,0), 
with n = 1.4 and «i = 1.28442 , trajectories are shown for 
the first 300 time units, the inner ring represents the orbit of 
electron 1 and the outer ring represents the orbit of electron 
2, the units are the scaled units defined in section II. 

FIG. 2. Non-ionizing double-ring tridimensional orbit for 
helium {Z = 2), trajectories are shown for the first 200 time 
units, the inner ring represents the plane projection of the 
orbit of electron 1 and the outer ring represents the projection 
of the orbit of electron 2. Positions are in the scaled units 
defined in section II. 

FIG. 3. Non-ionizing tridimensional orbit for H-minus 
(Z = 1), trajectories are shown for the first 10000 time units. 
The inner ring represents the plane projection of the orbit 
of electron 1 and the outer ring represents the projection of 
the orbit of electron 2. Trajectory of the (fastest) electron 1 
winds almost everywhere in the the dark inner core of figure. 
Positions are measured in the scaled units of section II. 

FIG. 4. Fast Fourier Transform of the orbit of Figure 2 
using 2 points. Frequencies are measured in the scaled units 
of section II. 

FIG. 5. Minimum time to ionization (among 24 random 
perturbations of average size 1 added to the orbit of Figure 
1) j3 is the adimensional parameter and time is measured in 
the scaled units of section II. 

FIG. 6. Average time to ionization (among 12 random per- 
turbations of average size 20/3 2 added to the orbit of Figure 3) 
(3 is the adimensional scale parameter and time is measured 
in the scaled units of section II. 
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